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Abstract. Simple modifications for higher-order 
Godunov-type difference schemes are presented which 
allow for accurate advection of multi-fluid flows in 
hydrodynamic simulations. The constraint that the sum 
of all mass fractions has to be equal to one in every 
computational zone throughout the simulation is fulfilled 
by renormalizing the mass fractions during the advection 
step. The proposed modification is appropriate for any 
difference scheme written in conservation form. Unlike 
other commonly used methods it does not violate the 
conservative character of the advection method. A new 
steepening mechanism, which is based on modification 
of interpolation profiles, is used to reduce numerical 
diffusion across composition discontinuities. Additional 
procedures are described, which are necessary to enforce 
monotonicity. Several numerical experiments are pre- 
sented which demonstrate the capability of our Consistent 
Multi-fluid Advection (CMA) method in case of smooth 
and discontinuous distributions of fluid phases and under 
different hydrodynamic conditions, ft is shown that 
due to the reduced diffusivity of the proposed scheme 
the abundance of some heavy elements obtained from 
hydrodynamic simulations of type 11 supernova explosions 
can change by a factor of a few in the most extreme cases. 

Key words: hydrodynamics - nuclear reactions, nucle- 
osynthesis, abundances - methods: numerical - super- 
novae: general 



1. Introduction 

One of the most important factors determining the quality 
of a numerical algorithm is its robustness. In the simplest 
case it can be regarded as a property of the scheme to pro- 
vide the result at minimum cost once the accuracy is spec- 
ified. In numerical hydrodynamics much effort has been 
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spent during the last two decades on improving advec- 
tion schemes in such a way that they do not only provide 
high accuracy in regions of smooth flows but also resolve 
flow discontinuities (shocks and contact discontinuities) 
with a minimum number of discrete grid zones. In the 
past a variety of numerical experiments were conducted 
to compare the overall quality of solutions obtained with 
the help of different advection schemes which provided an 
understanding of the advantages and disadvantages of al- 
ready avail able or newly propose d alg orithms (Col ella fc 
Woodward |1984 Carpe nter et al. |1990| , Fr yxell et al ~ 



Yang & Przekwas |1992| , Stone & Norman |1992| , Steinmctz 
& Muller 



1993, Kang et al. 1994) 



1991 



For our numerical experiments we have used the Piece- 
wise Parabolic Method (PPM) of Colella & Woodward 



1984; hereafter CW) to study the evolution of multi-fluid 



flows with strong discontinuities and stiff source terms. In 
theoretical astrophysics the PPM method has been used 
to study a range of hydrodynamic phenomena like stel- 
lar collisions (Ruffert & Muller |1990| , Frolov et al. |1994|), 



evolution of supersonic jets (Balsara & Norman 1992, Bas 



set & Woodward 1995), large-scale structure formation in 



cosmology (Bryan et al. 1994), interaction of stellar winds 
in massive close binaries (Stevens et al. 1992), and the 
stability of radiative shock waves (Strickland & Blondin 



1995). Numerical experiments of CW and Yang & Przek- 



was (1992) clearly demonstrated the superiority of the 



PPM scheme among several modern advection schemes. 

A particularly interesting and challenging astrophys- 
ical problem involving multi-fluid flow (and one of our 
numerical experiments; see section 3.4) is the simulation 
of mixing in supernova envelopes. Mixing occurs because 
the non-steady propagation of the shock wave formed af- 
ter core collapse gives rise to Ray leigh- Taylor instabili- 
ties (for a recent review, see e.g., Muller 1998). The first 
(2D) simulations of mixing involving ten separate fluids 
were performed by Arnett et al. (1989). They computed 



the propagation of the supernova shock wave, which was 
artificially created by a "point" explosion, through the 
envelope of a realistic stellar model. Better resolved and 
more detailed simulations were later performed by Fryxell 
et al. ( |1991| ) and Muller et al. ( |l99l[ ). They identified the 
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Rayleigh- Taylor unstable regions as being associated with 
discontinuities in the chemical composition in the envelope 
of the progenitor star. The simulations were performed 
with the PPM-based hydro dynamic code Prometheus, 
which keeps track of different nuclear species by solving 
a set of additional continuity equations (see Fryxell et al. 



1989) 



In the case of single-fluid advection the problem of dif- 
fusion across contact discontinuities plays a crucial role 
and provides a simple test case for studying mixing be- 
tween different fluids in numerical simulations. In multi- 
fluid flows both chemical and contact discontinuities may 
be present. Although, as we shall see later, there exist im- 
portant differences between both kinds of discontinuities, 
we nevertheless can profit from our experience of model- 
ing contact discontinuities when dealing with composition 
discon tinuit ies. In this context we point out that Fryxell 
et al. (|1989| ) demonstrated that the advection of a contact 
discontinuity is simulated better with a PPM scheme than 
with any other scheme they considered in their study. 

Mixing of different fluids cannot be ignored if the che- 
mical composition plays an important role in the hydrody- 
namic evolution. For example, in case of a realistic equa- 
tion of state the total gas pressure is calculated as the 
sum of partial pressures exerted by each kind of species. 
More complex physical processes, emission from an op- 
tically thin medium (radiative cooling) or absorption of 
radiation, are strongly sensitive to changes in chemical 
composition, especially to changes in the heavy element 
abundances. Last but not least, the process of nuclear 
burning, to which we will pay special attention later in 
this paper, directly depends on the amount and type of 
nuclear fuel. In this particular case, mixing of different nu- 
clear species due to numerical diffusion can substantially 
affect the final chemical composition or even the overall 
dynamics of the flow (for a recent review, sec e.g., Miiller 



1998 



The paper is organized as follows. In Sect, g we briefly 
describe the basic components of the PPM scheme and 
some specific features implemented into the Prometheus 
version used in our numerical experiments. We then give a 
detailed description of the new consistent multi-fluid ad- 
vection method. In section |^ we present the results of our 
test simulations. A discussion of the results is contained 
in Sect. 0. 



2. Numerical method 

In what follows, we consider only the Direct Eulerian 
formulation of the PPM scheme as implemented in the 



Prometheus code (Fryxell et al. 1989). Most of the pre- 
sented results are valid for and can efficiently be imple- 
mented in codes based on the Lagrangian with remap ap- 
proach, too. Both versions of the PPM method belong to 
the class of shock-capturing methods which are charac- 
terized by high-order spatial and temporal accuracy. In 



such schemes flow discontinuities are treated as solutions 
to the hydrodynamic equations (Riemann problem) and 
are represented by sharp (1 to 2 zones wide) and oscilla- 
tion free profiles of hydrodynamic variables. There is no 
need for adding large amounts of artificial viscosity for 
shocks to obtain correct post-shock states. Prometheus 



uses a Strang- type directional splitting (Strang 1968) for 
simulation of multidimensional flows. 

For given initial data and boundary conditions each 
hydrodynamic step of the PPM scheme begins with con- 
struction of the interpolants approximating the distribu- 
tion of flow variables inside each grid zone. The initial 
parabolic profiles are subsequently modified according to 
local flow conditions: density profiles are made steeper 
near contact discontinuities and the distributions of all 
variables are somewhat flattened near shocks in order to 
reduce high-frequency post-shock oscillations. Afterwards, 
monotonicity constraints are imposed on the interpolation 
profiles to avoid unphysical solutions. The monotonized 
profiles are used to calculate initial data for the Riemann 
problem at each zone interface by averaging the mono- 
tonized parabolae over the domain of dependence of the 
zone interface. The left and right states at the interface 
obtained in this way define the input data for the non- 
linear Riemann problem, which is solved iteratively. The 
solution of the Riemann problem provides average hydro- 
dynamic state variables at the zone interface, which are 
used to compute the fluxes for the advection step, whereby 
the hydrodynamic variables are advanced to the new time 
level. 

For simulations of mixing in supernova explosions the 
basic PPM algorithm was modified to allow for advec- 
tion of several nuclear species. The hydrodynamic state 
vector was extended by adding mass fractions for each of 
the species as were the left and right states used as in- 
put for the Riemann problem. The effective states at zone 
interfaces are obtained by averaging the mass fraction pro- 
files over a properly chosen domain of dependence for the 
zone interface. The interpolation step for the mass frac- 
tions also includes steepening and flattening procedures 
followed by a monotonization step. The mass fractions do 
not enter the Riemann problem. They are treated as pas- 
sive scalars and are advected with the flow depending on 
the upwind state as determined by the average velocity 
obtained from the solution of the Riemann problem. 

Prometheus includes the handling of a general equa- 
tion of state using the approach of Colella & Glaz (1985). 
Gravitational forces are included in the calculation of the 
effective states entering the Riemann problem and by a 
separate acceleration step at the end of each hydrody- 
namic sweep. Operator splitting is also used to couple nu- 
clear burning and hydrodynamics. The nuclear reaction 
network is solved using a multidimensional Newton iter- 
ation (Miiller 1986). A more detailed description of the 



imple menta tion of Prom etheu s can be found in Fryxell 
et al. (|1989|) and Miiller (|l998| ). 
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2.1. Multi-fluid advection 

We consider the one-dimensional initial-boundary value 
problem 

d t U + d x F(U) = G, 

with boundary data U{x — xl, t) and U(x — ir, t), where 
F, U and G are the flux vector, the state vector and the 
vector of source terms, respectively. In case of the Euler 
equations for multi-fluid ideal hydrodynamics the state 
vector is 



U = U(x,t) 




F(U) 



pu 
pu 2 + p 
(pE + p)u 
pX n u 



where p,u,E = e+it 2 /2, e,p, X n and pX n are the total gas 
density, the velocity, the specific total energy, the specific 
internal energy, the pressure and the mass fraction and the 
partial density of the n-th fluid, respectively. The closure 
relation for the above system is provided by the equation 
of state, p = p(p,e,X), where X is the vector of mass 
fractions. 

Let the number of different fluid phases be Nx and 
the mass fraction of the n-th fluid inside zone i be A^ jTl . 
Then the following relation must hold for t>to, 



(1) 



2.2. FMA method 



us- 



As it has been observed first by Fryxe ll et a l. ( 1989] ) 
ing PROMETHEUS and by Larrouturou ( |l991 ), the nonlin- 
ear character of higher-order Godunov-type schemes is the 
primary reason for the violation of (Q). During a simula- 
tion it is not guaranteed that in such schemes the sum of 
the mass fractions inside each zone remains equal to unity 
even if both the underlying advection scheme is conser- 
vative and the total mass of each fluid (summed over the 
whole grid) remains constant to within machine accuracy 
(see also Miiller |1998| ). 

This failure of high-order schemes can be understood 
when one realizes that interpolation profiles are con- 
structed independently for each mass fraction. Thus their 
sum can take an arbitrary value. We notice that this prob- 
lem has a predominantly local character: (i) it is most 
important in regions where changes in composition are 
substantial, and (ii) condition ([!]) can be violated for any 
subvolume of a zone (e.g. for the zone of dependence). 

In practice, condition (Q) is usually enforced by apply- 
ing a simple renormalization of the mass fractions after 
each step. This procedure, however, not only lacks any 
formal justification but it also leads to large systematic 
errors in the mass fractions of the least abundant species. 
It also violates the conservative character of the scheme. 



One possible solution to this problem is to modify the 
interpolation step in such a way that deviations of the sum 
of the mass fractions from unity will remain small inside 
each zone. According to the procedure proposed by Fryxell 
et al. (1989), which we will refer to as FMA, this can be 



obtained by calculating the sum of the mass fractions at 
the zone interface and to flatten the interpolation profiles 
for that zone totally, if 



A* - 



N x 

E 



X* - 1 



(2) 



with some predefined threshold value ££• Here, and in 
what follows, we take expressions involving # to mean # 
equal to either + or — , where + and — refer to the right 
and left interface of the z-th zone, respectively. 



In their original calculations Fryxell et al. (1989) used 
PROMETHEUS supplemented with FMA (with e E = 1CT 7 , 
while a less restrictive value of es = 1CP 3 seems to be ac- 
ceptable for most applications). Although FMA acts only 
locally, in the long term it affects large regions of the grid 
due to its high diffusivity which effectively reduces the 
spatial accuracy of advection of species to first order. Con- 
sequently, one has to expect a large amount of mixing of 
fluids especially near composition interfaces. 

2.3. CM A method 

The Consistent Multi-fluid Advection (CMA) method re- 
tains the high accuracy of the PPM advection scheme for 
mass fractions with constraint (Q) being accurately ful- 
filled, while any excessive flattening of interpolation pro- 
files (as seen in case of FMA) is avoided. 
The advection step can be written as, 

N 

y £ jPl X l {t + At)AV l = 

N N 

PiXimVi - At^2(AfF+ ATF-) (3) 

i=i i=i 

where Ff (with # 6 {+> ~ }) is the numerical flux vector 
across the zone interface of zone i. Af is the area of the 
zone interface, AVi is the zone volume and At the size of 
the time step. 

Since the PPM scheme is conservative, 

F7=fU, (4) 

the following condition holds for each component of pX 
separately: 



N 

E 

i=l 



PiXi^{t)AVi = const. 



(5) 



We seek for a set of numerical fluxes, J^f , which satisfy 
conditions (0) and (||) simultaneously. In general, the total 
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mass flux at a zone interface, pfuf, computed from the 
sum of the mass fluxes of individual species is not nec- 
essarily equal to the total mass flux computed with the 
(total) mass density pf: 



N X 



N x 



Y,F# n = '£pfufX* n ^pfuf. 

n—1 ri— 1 



(6) 



This inconsistency is the reason why any higher-order 
Godunov-type advection scheme in which the interpola- 
tion step for the mass fractions (or partial densities) is 
not appropriately constrained will violate condition (|l|). 
The inconsistency can be avoided when scaling the origi- 
nal partial mass fluxes in such a way that their sum 
is exactly equal to the total mass flux. 

We request the modified partial mass fluxes of zone i, 
J? (which depend on the modified mass fractions at the 

interface xf n ) to be consistent with the advection of the 
total mass, i.e. 



N X 



N x 



(7) 



Hence, we modify the original partial mass flux vector 
using the simple scaling operation 



(8) 



in such a way that it becomes consistent with the conti- 
nuity equation: 



ffJ2pf u f x t = pf u f- 

n=l 

Comparing (Q) with (Q) one obtains 



N x N x 

* V X* = V X* 
i / j A i,n / t ^i.m 

n—1 n—1 



(9) 



(10) 



i.e. the normalization constant, ipf, can be written in 
terms of the (unknown) modified mass fractions, Xf nl as 



A^/n—1 i,n 
t-^Nx x # 
Z—m— 1 i,7i 



(ii) 



Since the modified mass fractions are consistent with ad- 
vection of the total mass by definition (Q), their sum is 
equal to unity, i.e. 



1 



Z_-m— 1 i.n 



(12) 



In practice it is not necessary to explicitly compute the 
modified partial mass fluxes (ph, but instead one simply 
has to scale the original mass fractions X^ n according to 
©, i-e. 



X* = LD*X* . 

2,71 rz ^*~z,n 



(13) 



The CMA method does not require any modification of the 
interpolation step of the advection scheme. Furthermore, 
flux scaling does not destroy the conservative character 
of the scheme (^). Hence, scaling the average mass frac- 
tions ( |l3| ) obtained from the solution of the local Riemann 
problem by ipf defined in Eq. ( |T^ ) is sufficient to satisfy 
both (0) and (§. 

In passing we note that instead of scaling the partial 
mass fluxes one could, in principle, appropriately adjust 
the total mass flux, too (which is equivalent to comput- 
ing the total mass flux as the sum of the partial mass 
fluxes). Although numerical experiments (not presented 
here) demonstrate that this method produces results of 
comparable quality to the CMA method, we prefer to use 
the latter because it preserves the original role of the to- 
tal gas density in the PPM method. Finally, one could 
discard the continuity equation for the total density and 
use instead the partial densities as the primary variables 
to calculate the total density whenever needed. We found 
this method to give results of overall very poor quality, 
and being much more diffusive near contact and composi- 
tion discontinuities. 

Modification (|l^) implemented in the original PPM 
advection allows us to abandon the FMA scheme, and to 
retain the high accuracy of the PPM method. We will re- 
fer to this algorithm as sCMA. There are, however, still 
some problems left to be solved. In what follows, we fo- 
cus on the problem of numerical diffusion across compo- 
sition interfaces, and in these cases when the nonlincarity 
of the advection scheme becomes too strong to preserve 
the monotonicity of the scheme. This will complete the 
description of all of the elements constituting the CMA 
method. 

FMA uses a contact steepening mechanism which in 
PPM is used to limit diffusion across contact discontinu- 
ities only. Removing FMA and using the full capabilities of 
the PPM interpolation algorithm for mass fractions, how- 
ever, revealed a severe problem. Overshooting occurred 
near composition interfaces for the most abundant species. 
This, in turn, caused the least abundant species to be 
evacuated out of the critical region. Clearly, in its origi- 
nal version the PPM contact discontinuity detection algo- 
rithm cannot be directly used to steepen distributions of 
mass fractions. Some additional criteria have to be used 
to identify composition interfaces which need to be steep- 
ened and some mechanism has to be devised which limits 
overshooting. We point out that in most cases the FMA 
method effectively prevents any steepening of mass frac- 
tions profiles, since once any of the fluid distributions has 
been steepened, it likely activated the FMA flattening pro- 
cedure. 

Since we do not expect the geometrical properties of 
composition interfaces to be substantially different from 
that of contact discontinuities, most of the original steep- 
ening algorithm can be safely used without modification 
for detecting large gradients in the fluid's composition. It 
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is very unlikely to find some additional constraint similar 
to that originally proposed for detection of contact discon- 
tinuities (CW, Eq. 3.2) in order to make the scheme more 
selective. Therefore, we will only consider local proper- 
ties of composition profiles. One possibility is to associate 
composition interfaces with rapid changes in fluid compo- 
sition. For this purpose we define a steepness measure for 
mass fraction profiles 



X 



X 



i+2,n 



Additional steepening is applied if a^ n is larger than a 
fixed value ax- Since a^ n has a similar meaning as the 
parameter ui used in PPM for calculation of the flattening 
coefficients (CW, Eq. A8), we set ax = 0.75. To prevent 
the steepening of relatively small composition jumps the 
following criterion has to be satisfied, too: 



IX 



x. 



min(A"i 



Xi- Un ) > 0. 



We use e x — 0.01. Furthermore, we do not steepen com- 
position profiles near extrema. If the zone is located next 
to an extremum for which the following criterion has to 
be fulfilled 

(Xi + 2,n — Xi+l,n) (Xi-i in — Xj_2,n) < 0, 

the steepening coefficient is set to zero. Finally, we do not 
steepen abundance profiles inside contact discontinuities, 
because this gives rise to enhanced overshooting of the 
partial densities, piXi^ n . After steepening we flatten and 
monotonize the mass fraction profiles in the same way as 
in PPM. 

As a remedy for the overshooting near composition 
discontinuities, we introduce two additional modifications. 
Local extrema of the mass fraction distributions are iden- 
tified by the criterion 

{Xi+\,n ~~ X i n ) (X i n — Xj + x,n) < 0. 

If the criterion is fulfilled in a zone i, the interpolation 
profiles are flattened by an additional (constant) amount 
in the two zones next to it (i ± 1): 



X* = wX i±hn + (1 - w)Xf ±1 



# 



(14) 



We use w = 0.5. 

The second modification is somewhat more complex. 
For each zone we calculate two sums 

Sf' CT = £max(0,a(Xf n -X i> „)), o e {+,-}, (15) 

71=1 

at both zone interfaces (# 6 {+, — }), which give the total 
positive Sf' + and negative Sf'~ deviation between the 
values of the mass fractions at the respective zone interface 
xf n and the values of the zone averaged mass fractions 

Xi _ 77 . 



We then correct the interface values of the mass frac- 
tions according to (H). However, instead of the constant 

_LL 

flattening coefficient w, we use a variable coefficient wY n , 



which depends on the positive (Sf' 
deviations defined in (|l5|): 



and negative (S', t 



sf n max ( 0, min 11,/? 



A* 



A 



, (16) 



where 



min (Sf 



A 



# 

i.max 



= max(Sf' + ,Sf'~) 



and 



sign(X+ -X7 n ) #sign(sf+-^ 



with /3 = 0.25. 

The reasoning behind this procedure is based on the 
observation that any variation in the distribution of one 
fluid component should be compensated by an appropriate 
variation of the distribution of at least one other compo- 
nent. Here, we define two separate groups of fluid com- 
ponents which show deviations (from the zone average) 
of the same sign at the zone interface. Once the sums of 
the negative (S~) and the positive (S + ) deviations are 
obtained, we try to limit their relative difference. Note 
that the FMA flattening criterion is based on the abso- 
lute value of the deviation of the sum + S ,#,_ from 
unity. 

In the CMA method the amount of additional flat- 
tening is not constant but smoothly increases with the 
relative difference between the two deviations. Maximum 
flattening is applied if the relative difference between the 
two deviations exceeds 1//3, while no additional flatten- 
ing is introduced for /3 = 0. Hence, when (3 > 0, the in- 
terpolation profiles for that group of species which shows 
the largest absolute deviation are modified using the same 
flattening coefficient for all species. This reduces the de- 
viation and brings it closer to that of the other group of 
species. 

In this way we introduce a finite amount of coupling 
between the interpolation profiles of two distinct groups 
of fluid components rather than trying to adjust the dis- 
tributions of fluid components individually. Since the ad- 
ditional flattening procedure has a strictly dissipative cha- 
racter and relies on flattening of interpolation profiles, it 
is very unlikely that it causes unphysical solutions. 

Finally, we mention that the just described procedure 
could also be used to guarantee consistency between the 
total mass flux and the sum of the partial mass fluxes. 
Hence, it could be applied both in the interpolation and 
the advection step. 
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Fig. 1. Initial distribution of fluid phases in Sod's shock 
tube problem: X\ - solid line; X 2 - dashed line; X3 - 
dash-dotted line. 



Fig. 3. Long-term behaviour of mean and extreme de- 
viations from condition (Q) in Sod's shock tube problem. 
Top: original PPM; middle: FMA; bottom: sCMA. 



Fig. 2. Distributions of fluid phases (left) and deviations 
of the sum of mass fractions from unity (right) in Sod's 
shock tube problem at time t — 1 for three different runs: 
original PPM code (top), PPM with FMA (middle), and 
PPM with sCMA (bottom). X\ - solid line; X 2 - dashed 
line; X 3 - dash-dotted line. 

3. Results 

We have performed several numerical tests to illustrate 
problems arising in multi-fluid flows (when simulated with 
Prometheus) and to demonstrate the capability of the 
CMA method. 

3.1. Shock tube 

The first problem i s the shock tube test problem originally 
proposed by Sod ( 1978[ ). We have modified this problem 
to include three passively advected fluids. The initial state 
for this problem (Fig. |l|) is, 

U(0<x<0.5,t = 0) 

and 

[7(0.5 < x<1.0,i= 0) = I u 
with a discontinuous distribution of X\ . 

X 1 (x,t = 0) 




0.8, 0<a;<0.5, 
0.3, 0.5<x<0.75, 
0.1, otherwise, 



and oscillating mass fractions of the two other fluids, 

X 2 = 0.3sin 2 (207rx), 
X3 = 1 — X2 — X3. 

The simulations have been performed with an ideal gas 
equation of state with 7 = 1.4, and an equidistant grid of 
100 zones. Reflecting boundary conditions were imposed 
on both ends (x = and x = 1) of the computational 
domain. 

We obtained three sets of results for different code con- 
figurations: original PPM code (Fig. @, top row), PPM 
with FMA (middle section, e s = 10"^), and PPM with 
CMA (bottom part) but with no special modifications 
of the interpolation algorithm included (sCMA). The left 
panels of Fig. || show the distributions of the fluids {Xi 
- solid line, X 2 - dashed, X3 - dash-dotted). The right 



panels give the deviations of the sum of the mass frac- 
tions from unity, As (see Eq. (|])), with the lower right 
plot in Fig. U illustrating the truncation error of the CMA 
method (Eqs. @ and |[|)). 

The comparison of mass fraction profiles of different 
models at t = 1 (Fig. ||) shows a good agreement for 
x<0.85 except for some small amount of clipping near 
the extrema of X 2 and X 3 . However, we find large dif- 
ferences near the jump in X\ at x«0.87. The amplitude 
of the sinusoidal variation of X 2 and X3 is significantly 
reduced in case of FMA, and the initial discontinuities in 
X\ are strongly smeared. The original PPM scheme and 
the sCMA scheme are much less diffusive. Both jumps in 
X\ are clearly separated and the profiles of X 2 and X3 do 
smoothly vary near xpsO.87. 

The differences are even larger near the right bound- 
ary of the computational domain (Fig. ||). PPM strongly 
violates condition (Q). Yl-^i deviates from unity by up 
to 6% (where X\ is discontinuous). The errors are much 
smaller for FMA the deviations only being of the order of 
10 -4 which is close to the chosen value of £a- However, 
as already noted above, the smaller error is bought at the 
cost of a degraded resolution. The sCMA method gave 
the best result. It does not only advect the fluids with 
high accuracy, but it is also able to keep Yl-^i = 1 at 
the level of machine accuracy. The only imperfectness one 
notices is some overshoot in the distribution of X\ just to 
the left of the larger discontinuity signaling the first sign 
of a need for additional modifications of the interpolation 
scheme of the mass fractions (especially near composition 
interfaces) . 

To observe the long term behaviour of the three 
schemes we continued our simulations up to t ~ 400 (more 
than 75 000 steps with a Courant number of 0.8). Figure 
^ shows the evolution of the maximum negative and pos- 
itive deviations from condition (|l|) recorded for each time 
step together with the mean absolute value of the devia- 
tion from unity averaged over all zones. The results for the 
original PPM method (top panel in Fig. ||) indicate large 
variations (in excess of 20%) which would certainly destroy 
any solution sensitive to chemical composition. When us- 
ing FMA we observe a rapid rise of the maximum error 
which levels off after slightly exceeding 2 ea- The results 
obtained with sCMA show a slow growth of the maximum 
and minimum deviations, which seem to saturate at later 
times. 



3.2. Two interacting blast waves 

This t est problem was originally proposed by Woodward 
( |1982|) . It gained much of its popularity later when being 
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Fig. 4. Initial distributions of fluid phases in the interact- 
ing blast waves test problem: X\ - solid line; X2 - dashed 
line; X3 - dash-dotted line. 



Fig. 6. Long-term behaviour of mean and extreme devia- 
tions from condition ([!]) in the colliding blast waves prob- 
lem. Top: original PPM; middle: FMA; bottom: sCMA. 



Fig. 5. Distributions of fluid phases (left) and deviations 
of the sum of mass fractions from unity (right) in the 
interacting blast wave problem at time t = 0.038 for three 
different runs: original PPM code (top), PPM with FMA 
(middle), and PPM with sCMA (bottom). X x - solid line; 
X 2 - dashed line; X 3 - dash-dotted line. 



used by Woodward & Colella (|198J) in their study of vari- 
ous advection schemes in case of flows with strong shocks. 
In this problem the initial state consists of a low-pressure 
region located in the central part of the grid, 



C/(0.1<x<0.9,t = 0) 



which is bounded by two regions of (different) high pres- 
sure 





( 1 







; - 


I 0.01 








1 1 









V 




I 100 



U(0<x<0.1,t= 0) 



and 



U(0.9<x<l,t= 0) 



For our test runs we used three passively advected fluids 
with mass fractions that are initially varying smoothly 
across the entire grid (Fig. |l|), 

X x = 0.5a; 2 , 

X 2 = 0.5 sin 2 (20a;), 

X3 = 1 — X2 — X3. 

Again, we use an ideal gas equation of state with 7 = 
1.4. The grid consists of 400 equidistant zones. Reflecting 
conditions are imposed at both grid boundaries. 

The results of our simulations at t = 0.038 are shown 
in Fig. ||. 

Since this time the initial distributions of the fluids 
are smooth we do not expect to see any discontinuities 
in the distributions of the mass fractions at later times. 
However, a discontinuity is created in the Xi and X3 dis- 
tributions at x w 0.6 when using the original PPM and 
the sCMA method (left column, top and bottom panel of 
Fig. ||, respectively). No such discontinuity is present in 
the FMA data (middle panel). We identify the creation of 
such spurious composition interfaces with the actual fail- 
ure of the unmodified discontinuity detection procedure of 
PPM (see Section [Q| ) . FMA once again proves to be more 
diffusive than the other two schemes: high-amplitude vari- 
ations of X2 and X3 as seen in the PPM and sCMA data 



Fig. 7. Initial distributions of fluid phases for the shock- 
contact interaction problem: X\ - solid line; X2 - dashed 
line; X3 - dash-dotted line. 

for 0.65sxSi0.8 have markedly smaller amplitudes when 
calculated with FMA. Moreover, there is no trace of an 
extremum in X3 at x«0.75. In other parts of the grid all 
three methods produce very similar results. 

As in the case of Sod's shock tube problem the condi- 
tion (0) is most strongly violated when the original PPM 
method is used (upper left panel in Fig. ^). The maximum 
deviation of 2% occurs in that region where the FMA 
results are mostly affected by the use of an additional 
flattening procedure. On the other hand, FMA violates 
condition (|j) at the level of £a with a single pronounced 
maximum at the spurious composition interface created in 
the other two schemes. The sCMA method produces the 
most accurate results both during the initial phases of the 
evolution and in the long term evolution (lower left panel 
111 Fig. |). The maximum deviations from (Q) exceed 10% 
for the original PPM method and fluctuate between 2 and 
3 times ea in case of FMA. 

3.3. Shock-contact interaction 

The initial state for this problem is, 



U(0<x<0.1,t = 0) 





U(0.l<x<0.5,t = 0) 



and 



J7(0.5<x<l,i = 0) 



with Xi = 0.2 for x < 0.15, X x = 0.6 for x > 0.15, and 

X 2 = 0.3 sin 2 (30a;) + 0.001, 
X3 = 1 — Xi — X2 ■ 

The initial abundances with a composition interface be- 
tween Xi and X3 at x = 0.15 are shown in Fig. ^. Again an 
ideal gas equation of state with 7 = 1.4 is used. The grid 
contains 400 equidistant zones. The left grid boundary is 
reflecting, while a flow-in boundary condition is imposed 
at the right grid boundary. The state of the inflowing gas 
is equal to that of the gas located near that boundary at 
the initial time. 
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Fig. 8. Distributions of fluid phases (left) and deviations 
of the sum of mass fractions from unity (right) in the 
shock-discontinuity interaction problem at time t — 0.045 
for three different runs: original PPM code (top), PPM 
with FMA (middle), and PPM with sCMA (bottom). X x 
- solid line; X2 - dashed line; X% - dash-dotted line. 



Fig. 10. Density, velocity, temperature and composition 
profiles for the most abundant species (from top to bot- 
tom) of the 15 Mq progenitor model. The high tempera- 
ture region in the innermost (r<310 8 cm) part of the Si 
shell results from the deposition of 10 51 ergs of internal 
energy in order to initiate the explosion. 



Fig. 9. Long-term behaviour of mean and extreme de- 
viations from condition (^) in the shock-contact interac- 
tion problem. Top: original PPM; middle: FMA; bottom: 
sCMA. 

The initial conditions create a strong shock wave at 
x = 0.1 which propagates towards the left, hits the com- 
position interface (initially located at x = 0.15) and then 
collides with the strong contact discontinuity (initially lo- 
cated at x = 0.5) that slowly moves to the left. Upon 
interaction a pair of shocks is generated. 

Figure || shows the distribution of the mass fractions 
together with the deviations from the condition ([j]) at 
t = 0.045 after all strong interactions have already taken 
place. All three methods give comparable results in re- 
gions of pure advective flow (x <; 0.4). Towards the left 
follows a region with low-amplitude variations in the dis- 
tributions of X2 and X3, which is much more diffused 
when calculated with FMA (middle left panel in Fig. ||). 
The composition interface (x « 0.25) also seems to be 
smeared out in case of FMA, but it remains sharp in the 
other two cases. Finally, there is some overshoot in the dis- 
tribution of X\ which can be seen just to the right of the 
composition interface in the sCMA data. The interface is 
slightly broader when calculated with the sCMA method 
than with the original PPM method primarily because of 
the smooth rounded profile associated with the overshoot. 

The analysis of deviations from condition (Q) (right 
column in Fig. ||) confirms our conclusions from the pre- 
vious tests. Original PPM produces deviations of the order 
of a few percent. Most of the extra diffusion used by FMA 
to keep deviations small occurs in the region of interac- 
tion between the shock and the contact discontinuity, and 
near the composition interface. It is in this region where 
differences in abundances between FMA and the other two 
codes are most apparent. The long term behaviour of the 
deviations (Fig. ||) indicates that there is no tendency for 
deviations to become smaller with time. Using the original 
version of PPM we find errors in the several percent range 
(up to 10%). The typical deviations for FMA are between 
2 and 3 times e\. There is some systematic increase in 
deviations visible in case of sCMA at late times, but it 
does not seem to be of any importance. 



losion 



3.4- Supernova expl 

For our final "numerical" example we have chosen the 
shock propagation during the post-bounce evolution of 



a core collapse supernova. We consider the very early 
(i<3 s) stages of the evolution, when the just born super- 
nova shock begins to sweep through the stellar layers just 
outside the iron core triggering nuclear synthesis. We will 
focus our attention on the role which numerical diffusion 
plays in the process of nuclear burning and on its impact 
on the final chemical composition of the thermonuclear 
processed material. 

For these calculations we have used Prometheus in 
a version which allowed us to consider those physical 
processes which play an important role during the early 
phases of the shock propagation. The gas is assumed 
to be a mixture of 14 nuclei ( 1 H, 4 He, 12 C, 16 0, 20 Ne, 
24 Mg, 28 Si, 32 S, 36 Ar, 40 Ca, 44 Ti, 48 Cr, 52 Fe, 56 Ni). An 
a-network with 27 reactions coupling 13 of these nu- 
clei (all except 1 H) is used to describe nuclear burn- 
ing. The approximate equation of state includes contri- 
butions from radiation, the 14 fully ionized Boltzmann 
gases, and positron-electron pairs (included in an approx- 
imate way). Self-gravity of the stellar envelope is taken 
into account as well as the gravitational attraction of a 
central point source which mimics the nascent neutron 
star (M ns 1.28 M Q ). 

The simulations have been performed in one spatial di- 
mension assuming spherical symmetry. At the inner edge 
of the grid a reflecting boundary condition is imposed, 
while free outflow is allowed through the outer boundary. 
The inner boundary is situated at r in — 1.376 10 8 cm cor- 
responding to a mass coordinate of 1.317Mq, which is well 
inside the Si shell. The computational domain, which ex- 
tends out to a radius of 6.4 10 9 cm, is covered by 1600 zones 
(corresponding to a resolution of 40 km) in our standard 
resolution run. Additional simulations with up to 12 800 
zones (corresponding to a resolution of 5 km) have been 
performed to study the convergence behaviour of different 
advection schemes. 

The initial model is a 15 Mq pre-collapse model of 
a blue supergiant (S. Woosley, private communication), 
which closely resembles a progenitor model of SN 1987A 



(Woosley, Pinto & Ensman 1988). The density, velocity 



and temperature profiles of the initial model are shown in 
Fig. 0. 

In order to launch the supernova shock wave, we cre- 
ated a thermal bomb by adding E t b = 10 51 ergs in form of 
internal energy to the innermost (r < 3 10 8 cm) parts of 
the Si shell. The resulting shock wave propagates rapidly 
outwards heating up the stellar matter. Simultaneously 
a much weaker reverse shock propagates inwards. The 
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Fig. 11. Density (top), velocity (middle) and temperature 
(bottom) profiles of the progenitor at t = 3 s obtained with 
the CMA method and a grid resolution of Ar = 10 km. 
The supernova shock has reached the helium shell and 
hence has already passed several composition interfaces. 
At every interface a weak reflected shock is created, which 
can be recognized in the velocity profile. The large density 
jump near r ~ 10 9 ' 5 cm separates shocked gas that ini- 
tially formed the thermal bomb from matter in the stellar 
envelope. 

Fig. 12. Chemical composition of the ejecta obtained 
with the FMA method as a function of mass coordinate 
at t = 0, 100, 300, and 500ms (from top to bottom). The 
grid resolution is 40 km. 

Fig. 13. Chemical composition of the ejecta obtained 
with the CMA method as a function of mass coordinate 
at t = 0, 100, 300, and 500ms (from top to bottom). The 
grid resolution is 40 km. 



temperatures and densities behind the supernova shock 
are sufficiently high to trigger thermonuclear burning and 
the production of new elements . The resulting release 
of nuclear energy slightly enhances the explosion energy. 
Whenever the outward propagating shock crosses one of 
the composition interfaces of the progenitor star (see Fig. 
|l0| ), a weak reflected shock is created. The reflected shocks 
move inwards reheating the matter, which has just been 
processed by the supernova shock (see velocity profile in 
Fig. [ll]). A strong contact discontinuity at r w 10 9 5 cm 
corresponding to the mass coordinate M r w 1.415 Mq 
separates the shocked envelope gas from matter initially 
belonging to the thermal bomb. 



Figure 12 shows the evolution of the chemical compo- 
sition obtained with the FMA method at a resolution of 
40km in a narrow mass range (~ 0.1 Mq) close to the 
outer edge of the thermal bomb. Nuclear burning is most 
intense very early on (t < 100 ms) when matter is still 
dense and hot. At later times the chemical composition 
does not change appreciately with one exception. Produc- 
tion of 44 Ti only begins around t — 100 ms (a barely visible 
bump between 1.355 and 1.36 Mq) and lasts for ~300ms. 

In many aspects the evolution is similar when cal- 
culated with the CMA method at the same grid resolu- 
tion (Fig. |l^). However, some important differences also 
exist. With the CMA method mixing between 16 and 
other nuclear species at M r w 1.415 Mq is greatly re- 
duced. 36 Ar and 40 Ca are clearly separated from oxygen 
at t = 500 ms. There is an indication of a 28 Si interface 
near M r w 1.41 Mq. The transition region extends only 
over two zones (over about 5 zones in case of FMA) and is 
accompanied by a small amount of overshooting towards 
larger radii. This region is difficult to model due to the 



Fig. 14. Total mass of 40 Ca (upper part) and 44 Ti (lower 
part) as a function of time obtained with different advec- 
tion schemes. Solid lines: CMAZ (thin), FMA (medium), 
CMA (thick). CMA results with additional flattening for 
40 Ca are shown by dashed lines for fc a equal to 0.25 
(thin), 0.50 (thick), 0.75 (long thin), and 1.00 (long thick). 
The scale on the right side gives the masses normalized to 
the respective final mass obtained with CMA. 

Table 1. Total masses of 40 Ca and 44 Ti (in units of M ) 



at t = 3 s 


model 


40 Ca 


44 Ti 


CMA 


0.00718 


0.00050 


fo.25 


0.00699 


0.00067 


fo.50 


0.00684 


0.00085 


fo.75 


0.00670 


0.00105 


fl.OO 


0.00658 


0.00121 


FMA 


0.00662 


0.00122 


CMAZ 


0.00665 


0.00169 


WPE a 


0.006 


0.0002 



Model 15A of Woosley, Pinto & Ensman (198£ 



presence of the strong contact discontinuity separating the 
stellar envelope from the thermal bomb. Without addi- 
tional flattening and monotonization of the mass fraction 
profiles (as described in section ^) the overshooting of 
28 Si is suspiciously large. We note that overshooting of the 
most abundant species near composition interfaces might 
be a common problem for hydrodynamic codes (see, for 
example, Fig. 4 of Woosley, Pinto & Ensman 1988), and 
certainly deserves further investigation. 

Another difference between the FMA and CMA re- 
sults is the mass of 44 Ti produced in the simulations, 
which seems to be quite sensitive to the amount of nu- 
merical diffusion. Since titanium is synthesized via the re- 
action 40 Ca(a, 7) 44 Ti and since enhanced diffusion results 
in a smoother distribution of calcium, we computed sev- 
eral models where the interpolation profile for calcium was 
constructed assuming four different constant flattening co- 
efficients / Ca (models f . 25 , /0.50, Zo.75 and /1.00, respec- 
tively). An additional extreme model (CMAZ) was cal- 
culated where all mass fraction profiles are flattened com- 
pletely thereby imposing a maximum amount of numerical 
diffusion. The results (Fig. |l4|) show that the amount of 
44 Ti is smallest when using CMA (5.010~ 4 Mq), that it 
increases linearly with /c a and that the result of FMA 
(1.2210~ 3 Mq) is recovered when the calcium profile is 
kept totally flat throughout the simulation (model /i.oo)- 
In case of maximum numerical diffusion (model CMAZ) 
the amount of titanium (1.69 10 -3 Mq) is even larger and 
exceeds that obtained with CMA by a factor of more than 
three. Table summarizes these results; data taken from 
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Fig. 15. Composition structure of the progenitor at t = 
3s: CMAZ (top), FMA (middle), CMA (bottom). 

Fig. 16. Dependence of the production of 44 Ti on the grid 
resolution. The total 44 Ti mass is shown as a function of 
resolution at t — 3s for models CMAZ (open squares), 
FMA (open circles), and CMA (full circles), respectively. 



model 15A of Woosley, Pinto & Ensman ( 198§| ) (who used 
a different mechanism to initiate the explosion!) are also 
shown for comparison. 



Figure 15 shows the composition profiles in the ejecta 



at t = 3s for our three basic models: CMAZ (top), FMA 
(middle), and CMA (bottom). In CMAZ all abundances 
change smoothly and no particular feature can be recog- 
nized. In CMA there exist composition interfaces of 40 Ca 
(logr » 9.36), 16 and 32 S (logr rj 9.5), and several dis- 
continuities (in 28 Si, 32 S, 36 Ar, 40 Ca) at logr w 9.5. Out 
of these only a relatively weak discontinuity in the distri- 
bution of 36 Ar can be recognized in the FMA model. 

Finally, we have studied the convergence properties 
with respect to the production of heavy elements in our 
three schemes. The results for 44 Ti (Fig. |l6| ) indicate that 
with CMAZ and FMA the production of titanium de- 
creases as the resolution is improved. More interestingly, 
the CMA results depend only weakly on the resolution. 
This behaviour can be understood if we realize that there 
exists a composition interface in calcium near which ti- 
tanium is formed, and that (as demonstrated by our nu- 
merical experiments; Fig. pT[ ) the production of titanium 
grows with the diffusivity of the advection scheme. Once 
the composition interface is resolved and properly handled 
by the code, the final mass of titanium becomes practically 
independent of the spatial grid resolution. 



dynamic code Prometheus. Since the advection part of 
PPM is well tested for a single fluid, we have not consid- 
ered simple test problems with known solutions, like e.g., 
the linear advection of square profiles in the mass frac- 
tions. Instead we investigated the behaviour of the CMA 
method in case of ID test problems with both smooth and 
discontinuous composition profiles involving flows with 
strong hydrodynamic discontinuities. As for these prob- 
lems no analytical solutions are known, the correctness of 
the proposed method has been demonstrated by means 
of a convergence study. Although other methods converge 
to the same solution too, the CMA method is the only 
one, which simultaneously guarantees the mass constraint 
"^2 Xi = 1, i.e. the sum of the mass fractions is always and 
everywhere equal to one. 

In order to demonstrate the advantage of the CMA 
method, we have also studied a problem of astrophysi- 
cal relevance, shock-induced thermonuclear burning in a 
supernova explosion. It is shown that numerical diffusion 
near composition interfaces can change the composition 
of supernova ejecta by a factor of a few. The abundance 
of titanium is most severely affected in our test calcu- 
lations. The consequences of these findings for explosive 
nucleosynthesis calculations should be explored in more 
detail. 
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4. Summary and conclusions 

We have derived a numerical approach which guarantees 
that the sum of mass fractions equals unity in simulations 
of multi-fluid flows with higher-order Godunov-type meth- 
ods. Unlike other commonly used numerical methods, the 
proposed scheme preserves the conservative character of 
the underlying advection scheme. We would like to stress 
that even if the advection step is formally written in con- 
servation form, this does not necessarily imply that the 
scheme is conservative in case of multi-fluid flows. This 
fact is often overlooked. 

Modifications of the interpolation step are needed in 
higher-order Godunov-type methods to reduce the numer- 
ical diffusion near composition interfaces. These modifica- 
tions are described together with procedures that ensure 
the monotonicity of the scheme. 

The Consistent Multi-fluid Advection method (CMA) 
is proposed as a new method to accurately describe multi- 
fluid flows and is implemented in the PPM-based hydro- 
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